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Abstract 

Non-uniform grids and mesh adaptation have been a growing part of numerical simulation 
over the past years. It has been experimentally noted that mesh adaptation leads not only to 
locally improved solution but also to numerical stability of the underlying method. There have 
been though only few results on the mathematical analysis of these schemes (see for example 
Sfakianakis (2011)) due to the lack of proper tools that incorporate both the time evolution and 
the mesh adaptation step of the overall algorithm. 

In this paper we provide a method to perform the analysis of the mesh adaptation method, 
including both the mesh reconstruction and evolution of the solution. We moreover employ this 
method to extract sufficient conditions -on the adaptation of the mesh- that stabilize a numerical 
scheme in the sense of the entropy dissipation. 



1 Introduction 

Hyperbolic conservation laws appear in various applications. For example, fundamental physical 
laws, the conservation of mass momentum and energy, lead to the Euler equations of gas dynamics. 
Further examples arise in traffic flows, shallow water flows, magnetohydrodynamics and biology, see, 
e.g. LeVcquc (2002), Toro (2001), Godlewski & Raviart (1990), Freistuehler, Schmeiser & Sfakianakis 
(2011). 

Let us consider a scalar conservation law in one space dimension, 

Ut + /(«)«=0, x G [a,b], (6 [0,T], (1) 

with initial data uq £ L°°([a, b]). In order to simplify the presentation we assume e.g. periodic 
boundary conditions. 

Adaptivity is a main theme in modern scientific computing of complex physical phenomena. It 
is important to investigate the behaviour of adaptive schemes for hyperbolic problems, such as (1), 
which exhibit several interesting and not trivial characteristics. In this work we study the behaviour of 
certain geometrically driven adaptive algorithms when combined with the important class of entropy 
conservative schemes introduced by Tadmor (1987, 2003). 

In every time step t = t n the mesh we consider is: 

M™ = {a = xl < ■ ■ ■ < x% = b} 

with variable space step sizes hf — xf +1 — xf, i = 1 . . . N — 1; the mesh is reconstructed in every time 
step t n . Further, we consider a numerical approximation U n of the exact solution u over the mesh 
at time t = t n given as a 

cr = K,..., u n N }. 



The construction and evolution of our non-uniform meshes and the time evolution of the approxi- 
mate solutions is dictated by the Main Adaptive Scheme (MAS) which is described by the following 
procedure: 

• in every time step, construct new mesh according to the prescribed adaptivity criterion, 

• reconstruct the numerical solution over the new mesh, 

• evolve the numerical solution in time using the numerical scheme. 

MAS will be discussed in details in Section 2; we note here that the number of spatial nodes is fixed 
and that the reconstruction of the mesh is realized by moving its points according to the geometry of 
the numerical solution. 

The use of non-uniform adaptively redefined meshes, in the context of finite differences, was first 
studied, among others, by Harten & Hyman (1983), Dorfi & Drury (1987), and Tang & Tang (2003). 
The approach that we follow, for the mesh reconstruction step of MAS (Step 1) was first introduced by 
Arvanitis, Katsaounis & Makridakis (2001) and by Arvanitis (2002). Applications of MAS on several 
problems, point out a strong stabilisation property emanating from the mesh reconstruction Arvanitis, 
Makridakis & Tzavaras (2004), Arvanitis & Delis (2006), Arvanitis, Makridakis & Sfakianakis (2006), 
Sfakianakis (2011), Sfakianakis (2009). These stabilization properties led Arvanitis, Makridakis & 
Sfakianakis (2006) to combine MAS with the marginal class of entropy conservative schemes. The 
later were first introduced by Tadmor (1987) and further studied by LeFloch & Rohde (2000), Tadmor 
(2003), Lukacova-Medvidova & Tadmor (2009). They are semi-discrete numerical schemes which 
satisfy an exact entropy equality. On one hand these schemes are interesting on their own right, since 
they appear in the context of zero dispersion limits, complete integrable systems and computation of 
non-classical shocks. On the other hand they are important as building blocks for the construction of 
entropy stable schemes Tadmor (1987, 2003), Lukacova-Medvidova & Tadmor (2009). 

We note that classical techniques for the analysis of numerical schemes, only include the time 
evolution step of the procedure. In order though to have a complete picture of the quality of the 
numerical solution under a mesh adaptation procedure a broader analysis is needed. In this direction 
the work done in Sfakianakis (2009) has provided some constructive analysis tools. In the present 
paper though we are able to combine in one relation the effect of both the time evolution and the 
mesh adaptation. Let us point out that our approach allows us to represent the effects of both the 
adaptive mesh reconstruction as well as finite volume scheme in one conservative update relation over 
a reference uniform mesh (21). This approach allows to apply the entropy stability analysis and derive 
a sufficient mesh adaptation criterion to control entropy production. 



2 Main Adaptive Scheme (MAS) 

Using non-adaptive meshes -both uniform and non-uniform- the evolution of the numerical solution 
is dictated solely by the solution update. On the contrary, in the adaptive mesh case, two more 
phenomena need to be taken into account; the construction of the new mesh and the solution update. 
These steps comprise the Main Adaptive Scheme (MAS): 

Definition 2.1 (MAS). Given mesh M™ and approximations U n , 

1. (Mesh Reconstruction). Construct new mesh M™ +1 

2. (Solution Update). Reconstruct U n over M™ +1 to obtain U n . 

3. (Time Evolution). Evolve U n in time to compute U n+1 over Af™ +1 . 
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It is important to note that in the case of a fixed mesh, uniform or non-uniform, there is no need 
for mesh reconstruction (Step 1.) and effectively no need for Step 2. In such case MAS reduces to 
just the time evolution step (Step 3.) which is what is usually considered as a numerical scheme. 
The extra steps of the adaptive MAS, on one hand, change significantly both the computation and 
the analysis of the numerical approximations, and on the other hand are responsible for stabilization 
properties of the MAS. 

In the mesh reconstruction step (Step 1.) the mesh nodes are relocated according to the geometric 
information contained in the discrete numerical solution. The basic idea is geometric: 

in areas where the numerical solution is more smooth/flat the density of the nodes is low, in areas 
where the numerical solution is less smooth/flat the density of the nodes should be higher. 

In fact, the mesh reconstruction process can be chosen in any suitable way. One possibility is to 
use the monitor function which reflects the curvature of the numerical solution. For details we refer 
Sfakianakis (2009), Arvanitis & Delis (2006). 

Further, we consider the solution update procedure (Step 2. of MAS). The numerical solution U n 
is given on the old mesh M™ and is recomputed as U n on the new mesh M™ +1 . There are many ways 
for the reconstruction. In this work we use conservative piecewise constant reconstructions. 

Finally, for the time evolution step (Step 3. of MAS) we use any numerical scheme valid for non- 
uniform meshes. Denoting the mesh-solution pair by U n = {u™, i = 1 . . . N}, M" = {C™, |C"| = hf, i = 1 . . . N}, 
we obtain in the case of finite volume scheme 

< +1 = «? " )>p K-1/2) • (2) 

Here U n = {u™, i = 1 . . . N} is a reconstructed U n over M" +1 and the numerical flux F is decorated 
with " since it is computed over the updated values U. The numerical flux itself can be any numerical 
flux valid for non-uniform grids. We refer to Sfakianakis (2009) and Arvanitis & Delis (2006) for more 
details regarding both the implementation of numerical schemes over non-uniform meshes and their 
properties. 

2.1 Reference uniform mesh 

A schematic representation of MAS (2.1) in the form of mesh-solution pairs is the following 

{M™, U n ) mcsh adap S {M£+\ U n } num - schcmc ) {M x l+1 , U n+1 }, (3) 

where in the first part we have considered the Steps 1 and 2 of MAS and in the second part the Step 
3. 

In parallel to MAS and (3) we define a new set of mesh-solution pairs where the meshes are uniform, 
constant in time, of the same cardinallity as M™ and discretizing the same physical domain. 

Definition 2.2 (Reference uniform mesh-solution pair). Let {M X ,U} and {M, V} be two mesh- 
solution pairs with M x = {Ci,\Ci\ = hi, i = 1...N}, M = {C u \Ci\ = Ax, i = 1...N], U = 
{ui, i = i . . . N}, and V — {Vi, i — 1 . . . N}. We call {M , V} the reference uniform mesh-solution pair 
to {M X ,U} if 

• the meshes M x and M discretize the same physical domain, and 

• the following per-cell mass conservation is satisfied for every i = 1, . . . , N 

Axv i = h i u l . (4) 

We prove in Lemma 2.1 that the per-cell conservation property (4) is a result of a geometric 
conservation law. 
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Geometric conservation law 

Let us consider a a time dependent cell C(r) = (x±(t), X2(t)). We look for an appropriate conservation 
law, 

u t (x,t)+£(u(x,t),x) x =0, (5) 
that expresses the mass conservation of u over the moving cell C(r). Thus, by the Leibniz rule, 

— \ u(x,r)dx = u(x2(t),t)x' 2 (t) - u(xi(t) 1 t)x' 1 (t) + / u t {x,r)dx. (6) 

« T Jxi(r) Jii(t) 

If the mass of u over C(r) remains constant with respect to r, the following condition holds 

u t (x,T)dx = -u(x 2 (t),t)x' 2 (t) + u(xi(r),r)a/ x (r). (7) 

xi(t) 

Integrating (5) over C(r) we obtain 

Ut(x,r)dx+ / t)> x) x dx = 0. 

Ki(t) •'xi(r) 

Now, using (7) we get 

£(u(x 2 (t), t), x 2 (t)) - ^(u(a;i(r), r), xi(t)) = u(x 2 (r), t)x' 2 (t) - u{x 1 (t),t)x' 1 (t). 

A suitable flux function £ hence is 

£(u(x(t), t),x{t)) = u(x(t), t)x'(t). (8) 

Therefore, the strong formulation of (5) reads 

u t (x,r) + (u(x,T)x t ) x = 0, (9) 

which is referred in the literature as the Geometric Conservation Law (GCL), see e.g. Trullio & 
Trigger (1961), Thomas & Lombard (1979) , Cao et al. (2002). 

As previously announced we can attain the per-cell mass conservation property (4) by discretizing 
the corresponding GCL. 

Lemma 2.1. The per-cell mass conservation label (4) is a consequence of the geometric conservation 
law (9). 

Proof. For every given cell- value pair Cj, Ui and the respective reference pair Ci, Vi -as provided in the 
Definition 2.2- we set the moving cell C(t) = (xi(t), X2(t)) for t € [ri, t 2 ] to be a linear interpolation 
of Ci and Ci 

C(r) = T —^C % + ^—^Ci. 
t 2 - ri r 2 - ri 

Now, to attain a discrete version of (9) we integrate it over C(r) 

x 2 (t) pX2{r) 

u t (x,r)dx + / (u(x(T),T)x t (x,T)) x dx = 0, 

Ki(t) Jxi(t) 

and invoke (6) to get 



d_ 



/ u(x, r)dx — u(x 2 (t),t)x' 2 (t) +u(xi(t),t)x' 1 (t) + / (v,(x(r),r)xt(x,T)) x dx = 0. 

'xi(r) Jxi(t) 
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We discretize explicitly in [ti,T2], set At = t 2 — t\, and recall that \Ci\ = hi, and Cj = Ax to get 

1 , A , , £ 2 (t 2 ) - £2 (n) , xi(ra) - a;i(n) , £ 2 (t 2 ) - £ 2 (ti) xi(r 2 ) - xi(n) n 

— {Axv l -h i u l )-u i \-Ui Vui Ui = 

At At At At At 

or simply 

A:r Vi = hi Ui. 

□ 

Remark 2.1. The time variable r used in the previous proof refers to ficticious time; it does not 
correspond to the physically relevant time. 

Let us point out that in the theoretical analysis we will use the reference uniform mesh-solution 
pair to combine the effects of mesh adaptivity and the numerical update. 

In view of the Definition 2.2, and after applying the per-cell mass conservation (4) on the schematic 
representation (3) of MAS, we get for i = 1, . . . , N: 

h? < = Ax < mcsh adap S fc» + 1 = Ax fi? num - schcmc > h? +1 < +1 = Ax < +1 . (10) 

Now, by invoking (10) after multiplying with h" +1 , we can rewrite the scheme (2) over the uniform 
reference mesh 

At 

~Ax V'''+V2 ~ J, <-i/2, 
where F is the numerical flux function from (2) written in variables v. 



3 Entropy stablility 

Before stating the main theoretical result we introduce the following notations: 

A< +1/2 := < +1 ~ «?, (12a) 

B i+i/2 ■= . (12b) 

4 + 1/2 



.„ _ /«+i) + /K ")-2^; 1/2 

h+l/2 - Ay n 



(12c) 



i+1/2 

A ™+V2 ._ n+1 n {lOA} 

^ X i+l/2 — J 'i+l/2 x i+l/2'i l lzu J 

*W •= I ^ - J , (12e) 

Mf :=^(ffr_ 1/2 -F« +1/2 ), (12f) 
Now, we proceed with the main theoretical result. 

Theorem 3.1. We use the notations (12a)-(12f) and assume that the following condition holds 
At f/ „, At 



M * ~ 4A^i [ D *-V* " K " Ax {Rl - 1/2 + QU ' 2 + A - 1/2)2 ) ( A <-V2) 2 (I 3 ) 



+ (^"+1/2 ~ ^ |^(^+l/2 - QI+1/2 - A+l/2) 2 )(A«r +1/2 ) 2 |, 



(14) 
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where Ax, At are respectively the space and time steps that correspond to the numerical scheme (2). 
The mesh adaptation procedure (3) is used, where the corresponding reference uniform mesh is given 
in the Definition 2.2. Then the mesh adaptation procedure MAS (3) with the numerical scheme (2) 
for the time evolution step is entropy stable. 

Proof. The numerical scheme for the uniform variables reads, cf. (11) 



„,«+! _ f.n At ( pn 77';) 



We subtract v" to develop the respective incremental form 

At 



< +1 " v? = v? - < ~ 7^ (2^ +1/2 - 2F« 



i-1/2 



Equivalently 



At 
2A~a 



u i+i 



/(«?) - /K-r) + 2 ^- 1/2 « - <-r) | ( 1 - I 



( (X n +i /2 - Q? +1/2 ) A< +1/2 + (sti/2 + Q?-i/ 2 ) A<_ 1/2 ) . 



(15b) 



We point out that the term vf — vf is new and accounts for the mesh reconstruction and solution 
update steps of the MAS (2.1). 

We now express — vf in a conservative form with respect to Accordingly the size of Cj 

changes as: 

h- +i = k + A<+y 2 2 - A X tm (i6) 

and the the mass of u over Cj as: 

K +1 u? = h?i# - (A< + + 1 / 2 2 )_< + (A<£/ 2 2 ) + < +1 + (A/^V^! - (A<+y 2 2 ) + < (17) 
(17) recasts to: 

(Ax n+1/2 ) (Ax n+1/2 )^ (Ax n+1/2 ) (Ax n+1/2 )^ 

«*"«* jj« ««+" ^ »«-! + - ^ -.-I" ««■ 

This relation can also be written as a conservative difference 

(*<■£>- .,. ^-Py ^ A*C?>- : . CO . ^ (181 

= H i-l/2 - #1+1/2 ( 19 ) 

Replacing (19) in (15b) we obtain 

v i w i — U i-l/2 n i+\/1 



At 
2Ax 



( (i^+i/2 - Q?+i/ 2 ) A^Vi/2 + (b?-i/ 2 + Q?-i/a) - (2°) 
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which can be analogously written as a conservative update over the reference uniform mesh 

We note that the conservative difference i?jL]/ 2 — H™ +1 ^ 2 accounts for the mesh reconstruction and 
the solution update step of the MAS. 



In order to simplify the presentation of the rest of the proof we assume that the entropy and the 
conservative variables (v and v, respectively) coincide, i.e. we choose U(u) = \u 2 for the entropy 
function. Now, to recover the entropy-entropy flux representation of (20), we multiply it by the 
entropy variables vf, yielding 



where G is the numerical entropy flux. We have further following Tadmor (2003), 



(22) 



£ {x) = 



1 



A_ 1/2 Ai;?_i /2 + A+i/2Au? +1/2 



g (F E)< K* fAt\ 



< 



4 \Ax J 
K 3 ( At\ 2 



4 \Ax J 



(Bi+i/2 - Qi + i/2) 2 Aw 2 2 +1/2 + (Bi_ l/2 + Q z _ 1/2 ) 2 Avf_ 1/2 



(B i+ l/2 ~ Qi+1/2 ~ A+l/2) 2 Av i 2 +1/2 + (-Bj-1/2 + Q*-l/2 + A-l/2) & V i-l/2 



For more details see also Appendix A. Now, (22) reduces to 
At , „„ „„ \ . ,„ At 



_ UK) + — (Gr + i/2 " G?_ 1/2 ) < M? - ^ (A-i/ 2 A^_ 1/2 + A + i/ 2 A^ +1/2 ) 
(-Bi+1/2 - QI+1/2 - A+i/2) 2 Aw 2 +1/2 + (Bi_ 1/2 + Q*_ 1/2 + A-i/2) 2 Aw 2 _ 1/2 



if 3 (At 



+ 4 VAx 

Hence, the sufficiency condition for entropy stability reads 



MP < 



At 
AAx 



At 



A-1/2 - * A^A-1/2 + Qt-i/2 + A-1/2) 2 ) (Aw"_! 



At 



1/2; 



Remark 3.1. We point out that numerical scheme (21) 

„,n+l_.,n A* / Ax 1 cm ^ rr pn I 

^ - w >" ^ — I — U i+1/2 + r i+ i/ 2 - "TV #1-1/2 - b i-XI2 I ' 



Ax V At 



□ 



is written using uniform variables and incorporates both the time evolution step and the adaptation 
of the mesh. 

Example 1. To gain further insight into (16), (19) we refer to Figure 1 and provide three special 
examples. 
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l/" I 

"i-l 1 



I? 



it? i r? i i?. 



(a) The cell i moves to the left 



/)" 



(b) The cell i moves to the right 



< 




ay 








T 1 i, n 

1 1 '+ 1 


"i-i 1 8 






1 I'i-l" 




1 l i+l 1 ' i+1 1 


ft?-! 




h" 
"i+1 



(c) The cell i moves to both directions 



Figure 1: Three cases of cell movement. 



Cell moves to the left. 

This means that Ax"^ 2 = -rf , Aa; n +*/ 2 = -rf^, so = h\ - rf + rf_j = r?_ x + if. The 



L i+l/2 ' i ' "-"i-l/a 

mass of u satisfies h" +1 u" = rf^u^i + i"tt"; passing in v variables «" = Sr^i^li + 
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(23) 



• Cell moves to the right. 

This means that Ax^/j? = Ax?+}£ = /?, so = fc? + - if = rf + In this 



+ 1/2 — ^"j-1/2 — H ' DU "i — "i ~ H — ' i ~ 'i+1- 

case, the mass of u satisfies = rfu" + or in v variables vf = prt^ -r — — > , .... . ] - 

or 

m m 



i+i 



Cell moves to both directions. 

Similarly, Ax"^ 2 = = ~ r i-v 80 K +1 = K + + r l-i- Moreover h™ +1 i 
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rf-^ + hfuf + ll 



or v% 




i+l) 



or 




(25) 



Let us point that now the conditions (23)-(25) are simple enough in order to be tested for a moving 
mesh algorithm. We will report numerical experiments in a future work. 

4 Conclusions 

We provide in this work a new framework of studying the combined effect of mesh adaptation and 
time evolution of a numerical solution. This new method, has the benefit of being described over a 
uniform "underlying" grid that resolves the physical domain with the same number of discretization 
nodes as the numerical solution itself. To exhibit properties of this technique we study the dissipation 
of entropy due to the adaptation of the mesh. We have derived a sufficient condition for the mesh 
movement in order to guarantee that the overall procedure dissipates the entropy. Consequently the 
resulting numerical scheme MAS (2.1). 
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A Entropy stable schemes 

In this section, for the sake of completeness, we present parts -after some modifications- of the analysis 
conducted in Tadmor (1987, 2003). We moreover note that these works are to be viewed in the context 
of the seminal works on the subject by Friedrichs & Lax (1971), Mock (1978). 

A finite volume approximation of the one-dimensional conservation law is written in the form 

< +1 = < - H (*Wi/a(« n ) " *i-i/2«» • 



The corresponding viscosity form reads 

At 



u? +1 -u? = 



2Ac 



-(/« + i)-/«))+Qr + i/ 2 A< +1/2 

- (/«) - /«_!)) - Q?_ 1/2 Au?_ 1/2 ] (26) 



with the numerical viscosity given by Q" + i/2 = ^a«'° + ^ 2F "+ 1 / 2 _ After denoting B i+1 / 2 



^~Au" " e< l ua tion (26) can be rewritten as 



4+1/2 



2Ai(( B,+1/2 ~~ <9»+i/2)a< + i/ 2 + (-8,-1/2 +Qt-l/2)A<_ 1/2 



Now, setting D i+l / 2 = Qi+1/2 — Qi+1/2' with Q* the numerical viscosity of an entropy conservative 
scheme we obtain 



,71+1 n At 



"I" ' ~ - «<' = ((^+1/2 - QI+l/2 - A+l/2) A<+i/2 + (Si-l/2 + Ql-1/2 + A-1/2) A<-l/2 

Next we multiply the last equation with the entropy variables v. After noting that 

(«r,< +1 - o = t/« +1 ) - c/«) - fr £) (c +1/2 ), 

where ^^'(v"" 1 " 1 / 2 ), £f (u") are given by 

:-(^)f„-;™+l/2N _ Z^" 1 / 2 /l A / „+l/2 rr/-n+l/2w -n+1/2 



£=-V 2 



- J (A«r_ 1/2 , J D?_ 1/2 Ae i/2 ) + J (A5? +1/2 , D i+1/2 Av? +1/2 ) , 



we arrive to the entropy-entropy flux pair representation 



C/« +1 )-C/«) + |^ 



{^ +1/2 - ^-1/2} = Ef E \v^) 



Now, entropy stability implies that: 



ef E \v n+1 l 2 )-^Sf{v n )<Q. (27) 
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In order to satisfy (27) the following conditions are sought as sufficient 

K3 ^ x ^ Bi + 1 / 2 ~ Qi+l/2 ~ A+l/2) 2 < A+1/2 
K ^ {Bl - 1/2 + Q ^- 1 ' 2 + A -!/2) 2 < A-l/2 

for every i. Due to symmetry they give 

R3 ^C iBl+1/2 ± Q '* +1 / 2 ± A +V2) 2 ^ A+l/2 (28) 

The last relation yields the following results for c = J^^ t 

A+1/2 > 
(B±Q* ± A+1/2) 2 <cA+i/2- 

This leads to Z) 2 +1/2 + (2(Q*±B)-c)A+i/2 + (<3*±5) 2 < 0. By setting ki = Q*+B and k 2 = Q*-B 
we get two different quadratic inequalities that need to be satisfied for A+1/2: 

A+1/2 + (2^1 - c) A+1/2 + fc 2 < and Z? 2 +1/2 + (2fc 2 - c) A+1/2 + fc 2 2 < 0. 

The necessary restrictions for the existence of a solution are c > 4fci and c > 4/c2. This mean that 
Q* ±B < c/4, i.e c > 4Q* or 

At 1 
Ax ~ 4K 3 Q* ' 
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